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Abstract We present a new approach for the simulation of Beyond Standard Model (BSM) physics within 
the Herwig++ event generator. Our approach is more generic than previous methods with the aim of 
minimising the effort of implementing further new physics models. Spin correlations, which are important 
for BSM models due to new heavy fermions and bosons, are discussed and their effects demonstrated 
for the Minimal Supersymmetric Standard Model (MSSM) and Randall- Sundrum Model using our new 
framework. 
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1 Introduction 

In the arena of modern particle physics Monte Carlo event 
generators have become essential tools for analysing ex- 
perimental data. They are necessary in order to compare 
the behaviour of theoretical predictions under the con- 
ditions present within a collider experiment by giving a 
realisitic description of the final-state particles which in- 
teract with the detector including any experimental cuts. 
It is essential that these generators reproduce Standard 
Model physics as accurately as possible since these pro- 
cesses will provide a background to any new physics sig- 
nals that might be present at future colliders as well as 
being of interest in their own right. New physics mod- 
els will also need to be incorporated into a Monte Carlo 
simulation in order for their implications to be fully un- 
derstood. 

There are a wide variety of new physics models and 
while one could implement each model independently in 
its own event generator it is more efficient to have a gen- 
eral purpose event generator that can handle a variety of 
these models but can also offer the full event simulation 
framework, i.e. hard process, decay, parton shower and 
hadronization. This is the approach that will be described 
in this paper using the Herwig++ [1,2] event generator. 

Herwig++ is a new event generator, written in C++, 
based on the well tested HERWIG [3-5] program. It is 
not simply a translation of the old FORTRAN code in 
to CH — h, it includes significant improvements to both the 
physics models and simulation framework. The object ori- 
ented aspect of the C++ language will allow future ad- 
ditions and modifications to be incorporated more easily. 
One area where improvements are needed is in the simu- 
lation of Beyond Standard Model (BSM) physics. In the 



past each model was hard coded in the generator making 
the addition of new models a time consuming process. We 
wish to minimise the effort required in order to add new 
physics models to Herwig++. Our approach is to factor- 
ize the problem in to smaller pieces and reuse as much 
information that has already been calculated as possible. 

For example in gluino production our method would 
first calculate the 2 — > 2 production matrix element then 
choose a decay mode for the gluino based on the branch- 
ing ratios and generate the decay products. In addition, 
since the gluino is a coloured object, there will be QCD 
radiation which is simulated more easily in our factor- 
ized approach as it is simply another step between the 
production and decay of the particle. While factorizing 
the problem in this manner makes many things easier it 
docs introduce complications when considering spin cor- 
relations. Additional information must be passed between 
each step to ensure that the final decay products are cor- 
rectly distributed including correlations between the pro- 
duction and various decays. 

Other packages, such as MadGraph [6], CompHEP [7], 
Shcrpa [8] and Omega [9] with WHiZard [10] exist which 
are capable of producing a wide class of BSM physics pro- 
cessed but they have limitations. The main problem is the 
efficiency with which the variety of possible processes can 
be generated. The above programs all treat the processes 
as 2 — > n scattering which requires the exact production 
and decay chain to be specified from the beginning. While 
this does mean that effects such as spin correlations are 
included automatically, it limits the number of processes 
that can be generated in a reasonable amount of time. 



1 A general list of programs for BSM physics can be found 
at htt p : / / www .ippp.dur.ac.uk /mont ecarlo /BSM/ 
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For example in order to study two different decay modes 
of the gluino, with the above generators, one would have 
to calculate the production step twice whereas our method 
would simply be able to pick another decay mode without 
recalculating the hard subprocess. 

To minimise the amount of work needed for every new 
model, our approach will require only a set of Feynman 
rules, the specification of any new non-Standard Model 
particles and their properties. In addition there is a mech- 
anism to read in parameters from a Les Houches [110 
file for a supersymmetric (SUSY) model. The new physics 
models currently implemented are the Minimal Super- 
symmetric Standard Model (MSSM) with CP, R-parity 
and flavour conservation and a Randall-Sundrum [12] type 
model with the lowest Kaluza-Klcin excitation coupling to 
Standard Model matter. 

Section [2] introduces how spin correlations are dealt 
with since these will be important when dealing with heavy 
fcrmions and vector bosons. Some details on the techni- 
cal structure of the code will be given in section [3] along 
with some physics of the models implemented thus far. A 
comparison with the FORTRAN code for some physical 
distributions will be presented in section[4]to demonstrate 
the consistency of our approach. 



2 Spin Correlations 

Many new physics models predict the existence of new 
particles that are as yet undetected by experiment. Heavy 
spin- 1 , spin-1 and spin-2 particles will be produced which 
will decay to lighter states. Their non-zero spin gives cor- 
relations between the production and decay steps which 
must be taken into account in order for the final-state an- 
gular distributions to be correct. An algorithm for dealing 
with these correlations is demonstrated in [13-15]. It will 
briefly be described below for the process e + e~ — ► tt where 
the top quark subsequently decays, via a W-boson, to a b 
quark and a pair of light fermions. 

Initially the outgoing momenta are generated accord- 
ing to the usual cross-section integral, 
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is the matrix element for the initial hard 



process and A t j are the helicity of the t and t respectively. 
One of the outgoing particles is then picked at random, 
say the top and a spin density matrix calculated 
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with N defined such that Tr p = 1 . 

The top is decayed and the momenta of the decay prod- 
ucts distributed according to 

(2tt) 4 f d 3 p b d 3 p w + 
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where the inclusion of the spin density matrix ensures the 
correct correlation between the top decay products and 
the beam. 

A spin density matrix for the W + is calculated because 
the b is stable 



p w+ , 



bW + 

W + 



(4) 



and the W decayed in the same manner as the top with 
the inclusion of the spin density matrix here ensuring the 
correct correlations between the W decay products, the 
beam and the bottom quark. 

The decay products of the W are stable fermions so 
the decay chain terminates here and a decay matrix for 
the W 



bW 

W+ 



M 
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(5) 



is calculated. Moving back up the chain a decay matrix 
for the top quark is calculated using the decay matrix of 
the W 7 



A t A' 



Currently only SLHA1 is supported. 



A w+ A w+ ' 

Since the top came from the hard subprocess we must now 
deal with the i in a similar manner but instead of using 
5\ i \ i when calculating the initial spin density matrix, the 
decay matrix of the top is used and the i decay is generated 
accordingly. The density matrices pass information from 
one decay chain to the associated chain thereby preserving 
the correct correlations. 

The production and decay of the top, using the spin 
correlation algorithm, is demonstrated in figures[T][2]and[3l 
The hard subprocess and subsequent decays were gener- 
ated using our new method. The results from the full ma- 
trix element calculation are also included to show that 
the algorithm has been correctly implemented. The sep- 
arate plots illustrate the different stages of the algorithm 
at work. Figure [1] gives the angle between the beam and 
the outgoing lepton. The results from the simulation agree 
well with the full matrix element calculation which demon- 
strates the consistency of the algorithm for the decay of 
the t. 

Figure [2] gives the angle between the top quark and 
the produced lepton. This shows the same agreement as 
the previous figure and demonstrates the correct imple- 
mentation of the spin density matrix for the t decay. Fi- 
nally figure [3] gives the results for the angle between the 
final-state lepton/anti-lepton pair showing the correct im- 
plementation of the decay matrix that encodes the infor- 
mation about the t decay. Again there is good agreement 
between our numerical results and the full matrix element 
calculation giving us confidence about the implementation 
of the spin correlation algorithm. 

The above procedure is well suited for implementa- 
tion in an event generator, as demonstrated, where one 
would like additional processes to occur between the hard 
production and decay such as showering of a coloured par- 
ticle. The algorithm as presented here is implemented in 
HerwigH — h and will be used extensively during the simu- 
lation of many BSM physics models. 
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Figure 1. Angle between the beam and the outgoing lepton in e + e~ — » ti — » bbl + vil~i>i in the lab frame for a centre-of-mass 
energy of 500 GeV with (a) unpolarised incoming beams, (b) negatively polarised electrons and positively polarised positrons 
and (c) positively polarised electrons and negatively polarised electrons. 
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Figure 2. Angle between the lepton and the top quark in e + e~ — > ti — > bbl* v\l~ v\ in the lab frame for a centre-of-mass energy 
of 500 GeV with (a) unpolarised incoming beams, (b) negatively polarised electrons and positively polarised positrons and (c) 
positively polarised electrons and negatively polarised electrons. 



3 Technical Details 

Instead of following the paradigm of implementing a spe- 
cific model we have chosen a more generic approach to 
the problem, which is intended to be as model indepen- 
dent as possible. This will allow a wider variety of models 
to be implemented within the event generator framework. 
The structure of the code relies heavily on the inheritance 
facilities available in the C++ language which allow inde- 
pendent structures to have a common heritage. 

Due to the existing structure it was sufficient only to 
consider generation of the hard subprocess and decay of 
the subsequent unstable particles since the showering and 
hadronizing are handled (almost) independently of the 
model details. Both the hard process and decay require 
a knowledge of the Feynman rules, couplings and masses 



within the model and these are currently implemented for 
the MSSM and Randall-Sundrum model. 

In Herwig++ the Feynman rules are encoded in Ver- 
tex classes. They form part of the structure that enables 
the calculation of matrix elements using the HELAS for- 
malism [16]. These classes provide the couplings to the 
particles defined within the model and so they must be 
provided with every new model which is implemented. The 
way in which the vertices have been set up minimises this 
effort and will be described in the next section. 



3.1 Vertices 

For a given combination of spins interacting at a vertex, if 
we assume the perturbative form of the interaction, there 
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Figure 3. Angle between the outgoing lepton and anti-lepton in e + e~ — > ti — > bbl + vil~ i>i in the lab frame for a centre-of-mass 
energy of 500 GeV with (a) unpolarised incoming beams, (b) negatively polarised electrons and positively polarised positrons 
and (c) positively polarised electrons and negatively polarised electrons. 
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Figure 4. Feynman rule for the Xi X°jZ° vertex where the O t 
are a combination of neutralino mass matrices. 



is a specific Lorentz structure and a limit on the num- 
ber of possible couplings for any given interaction. This is 
carried into the implementation of the vertices by defining 
a base class that holds all common functionality and in- 
heriting from this class to define a specific spin structure. 
The spin structure must then be further specialised into 
the exact vertex required by specifying the coupling and 
the particles that are able to interact at it. 

As an example, consider the XiX'jZ vertex in the 
MSSM as shown in figured 

A more general rule than that given in the figure is: 



Y(1-7s) + -2-( 1 + T*) 



(7) 



where c is the overall normalisation and 

a {L,R} are the 

left and right couplings respectively. We would choose c = 
g/ cos9 w , a L = and = Oj*. The overall coupling 
c is stored in the base class and the inherited class for 
the specific spin stores the left and right couplings, since 
they may not always be required. Finally the actual vertex 
class implements a function to calculate the value of the 
couplings. 



In addition to storing the couplings, the spin specific 
vertices have functions that can be used to either evaluate 
the vertex as a complex number or return an appropriate 
off-shell wave function. The ability to calculate not just 
the entire vertex but off-shell components underlies the 
HELAS [16] formalism for the calculation of matrix ele- 
ments. As an example consider decay of the top as in sec- 
tion [2] The HELAS approach factorizes the problem into 
two parts. First a vector wave function for an off-shell W + 
is calculated for a specified helicity of the top and bottom 
quarks. This is used as an input, along with the spinors 
for the light fermions, at the second vertex to calculate 
the final matrix element for that helicity combination. To 
obtain the spin-summed matrix element the procedure is 
repeated for all possible helicities of the external particles. 
This dramatically reduces the amount of code required for 
numerical evaluation of the vertices. It is also has the ad- 
ditional benefit of providing basis stated for the particles 
that can be stored and passed between the production and 
decay to ensure the spin correlations are consistently im- 
plemented. Our implementation and the formulation are 
described in more detail in appendix [A"l 



3.2 Decayers 

The decay of the top described in section [2] is handled by 
a class that is solely responsible for this decay. When just 
considering the Standard Model this is reasonable since 
the top is the only heavy fermiorQ present. However, in 
the case of new physics models there will be a wealth of 
heavy particles which decay and creating a class for each 
decay would be inefficient. Instead we have implemented 
a set of classes for each set of possible external spin states 



3 In this case the spinors for the quarks and light fermions. 

4 Heavy in this context means that it decays before it 



hadronizes. 
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and each of these classes, called Decayers, is responsible 
for a specific external spin configuration in a decay. At the 
present time only two body decays have been considered^. 
Appendix [C] contains a list of the currently implemented 
Decayer classes. 

The standard way in which Herwig++ handles particle 
decay modes is with a text file listing each mode, along 
with branching ratios and the object that will handle the 
decay. We have instead taken the approach of constructing 
the decay modes and Decayer objects automatically. All 
that is left to the user is to specify the particle(s) which 
will be decayed. The steps for creating the decay modes 
and Decayer objects are: 

1. specify the particles for which decay modes are re- 
quired; 

2. analyse each vertex to find whether the particle can 
interact, if it cannot skip to the next vertex in the list 
until one is found that is able to; 

3. find the decay products and test whether the decay 
would be kincmatically possible, if not skip to the next 
possible mode; 

4. if an object already exists that can handle the decay 
then assign it to handle the mode else create a new 
Decayer and assign this to handle the decay. 

The created Decayer object contains the appropriate code 
for calculating the matrix element for all possible helicity 
combinations which can be used in the spin correlation 
algorithm from section [2l 

There is an additional point to consider when dealing 
with SUSY models. These models contain additional pa- 
rameters, such as mixing matrices, that are necessary for 
vertex calculations. A mechanism has been implemented 
to read this information from a SUSY Les Houches Accord 
file [11] . In principle the file can contain decay modes along 
with branching ratios. If this is the case then the decay 
mode is not created automatically it is just assigned an 
appropriate Decayer. 

In order to be able to decay the particles we must first 
produce them along with their associated momenta. The 
next section will describe how this procedure is accom- 
plished within the new framework. 

3.3 Hard Processes 

In the hard process the initial momenta of the outgo- 
ing particles from a hard collision are calculated via a 
leading-order matrix element as in equation [T] There will 
be additional PDFs involved if the incoming particles are 
composite. In HerwigH — h the mechanism for this is again 
"factorized" into pieces concerned with the phase-space 
evaluation and pieces concerned with the calculation of 
|A4| 2 . Here we only concern ourselves with the 2^2 
cross sections and as a result the existing structure only 
requires us to implement calculations of \A4\ 2 for the new 
processes. 

5 Any three body decays read from a decay table are handled 
with a phase-space decayer and therefore will not include spin 
correlation informatation. 



Figure [5] shows the possible topologies for a 2 — > 2 
process at tree level. Obviously some processes will not 
involve all topologies. As in the case for the Decayers our 
approach is not to create a class for each possible process 
but instead create classes based on the external spins of 
the particles involved. The user simply specifies the in- 
coming states and an outgoing particle, all of the possi- 
ble diagrams with this outgoing particle are then created 
along with the appropriate MatrixElement object. The ob- 
ject is responsible for calculating the spin-averaged matrix 
element. In the HELAS approach the matrix element is 
first calculated by computing the complex amplitude for 
each diagram of a given helicity combination. The diagram 
contributions are then summed and the modulus squared 
taken. This is done for each helicity and the sum of each 
helicity gives the spin-summed \M\ 2 . In the case of strong 
processes it is easier to separate the colour structure from 
the evaluation of the matrix element. In a given process 
the diagrams can be split into "colour flows" , which are a 
combination of diagrams with the same colour structure, 
reducing the amount of computation required. 

To demonstrate this procedure consider the process 
gg — ► gg. The diagrams that contribute are shown in fig- 
ure[6l If the amplitude, stripped of the colour information, 
for the ith diagram is denoted by M.i the full amplitude 
for the 3 diagrams is given by, 



9i = 


if aci if Hd Mi, 


(8a) 


92 = 


if adi if blc M 2l 


(8b) 


93 = 


-if aib if icd M 3 , 


(8c) 



where gi denotes the full amplitude and / denotes the 
anti-symmetric structure constants. 

The factors of % associated with the structure constants 
are present because of the way the vertex rules are defined 
within the code. The vertices in Herwig++ are stripped of 
their colour information and therefore require the extra 
factors of i to be included, where appropriate, with the 
colour matrices to give the correct sign. 

The combination of structure constants in equation [8cl 
can be rewritten using the Jacobi identity to give 

.93 = (f aci f Ud - M 3 , (9) 

making it apparent that the colour structure of the s- 
channel gluon exchange diagram is simply a combination 
of the other two colour structures. The full colour ampli- 
tude can therefore be written as 

Mt = -[ci(M 1 -M 3 )+c 2 (M 2 + M 3 )}, (10) 

where Cj denotes the combination of structure constants 
from above and the combination of diagram amplitudes 
will be known as "colour flows" denoted by /j. Note that 
the overall minus sign can be dropped since it simply cor- 
responds to a phase that will not contribute to the final 
answer. In order to calculate \Ai\ 2 the constants Ci need 
to be squared. For the process being considered \c\\ 2 = 
\c 2 \ 2 = N 2 (N 2 -l) andcic| = c 2 c\ = N 2 (N 2 -l)/2 where 
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Figure 5. Tree-level topologies for a 2 
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2 process. The arrows denote the flow of momenta. 
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Figure 6. Diagrams contributing to the process gg — > 
where lowered Greek letters denote space-time indices and 
raised Roman letters denote colour indices in the adjoint rep- 
resentation, u, v are the spinors for the gluinos and e M) „ are 
the polarisation vectors of the gluons. The momenta (pi,P2) 
are incoming and (pz,pi) are outgoing. 



N c is the number of colours. The spin-summed matrix el- 
ement, averaged over initial colours and polarisations, is 
then 



1 1 



1^1 =o7 



1 



2 4(AT2 



o,/;7; A 



(ii) 



where Cij is a matrix containing the squared colour factors 
and /j A is the ith colour flow for the set of helicities A. 

As well as calculating \Ai\ 2 for a given process each 
MatrixElement object is also responsible for setting up the 
colour structure of the hard process which is required to 
generate the subsequent QCD radiation and hadroniza- 
tion. Depending on the colours of the internal and ex- 
ternal states involved there may be more than one pos- 
sible colour structure for each diagram. For the example 
in figure [6] each diagram has 4 possible colour topologies 
since the both the gluon and gluino carry a colour and 
an anti-colour line in the large N c limit. When an event 
is generated, if this case presents itself, a colour structure 
is picked at random from the N t possiblities. While this 
does rely on using the large N c limit we believe it to be 
an good approximation for the colour structures that we 
are dealing with. 

The possible Majorana nature of external states also 
gives rise to complications when calculating the matrix 
element. If the incoming states are a spinor and a barred- 
spinor then in the case where a u-channel diagram is re- 
quired two additional spinors must be calculated. The rea- 
son for this is that, using the notation of figure [5l when c 
and d are crossed their fermion flow can no longer be re- 
versed since the initial fermions set the direction in which 
these arrows point. The two additional spinors required 
are a spinor for the original outgoing barred state and a 
barred spinor for the original outgoing spinor state with 
care being taken to associate the new spinors with the 
correct helicity. 
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Figure 7. Resonant graviton exchange from gluon fusion to 
produce 2 fermions. 

Appendix [B] contains a list of the currently imple- 
mented MatrixElement classetQ. The example given above 
demonstrates a possible SUSY gaugino production process 
that is taken into account with our new mechanism. We 
have also implemented another mechanism for the simu- 
lation of resonances. This will be described now using an 
example from the Randall- Sundrum model. 

3.4 Resonant Processes 

Often we are interested in the study of s-channel reso- 
nances which decay to Standard Model particles rather 
than the production of a new particle in a 2 — > 2 pro- 
cess. We therefore include a mechanism to study this type 
of proccess. We will take as an example here the virtual 
exchange of a graviton, the lowest lying state of a Kaluza- 
Klein tower. The graviton is predicted by various models 
with extra dimensions where gravity is allowed to propa- 
gate in the bulk, an example process is shown in figure [71 
The same matrix element classes which are used to calcu- 
late the hard processes are used to calculate the resonant 
processes. Now, however, there is less computation since 
each MatrixElement will contain only a single s-channel 
diagram and hence a single colour flow. 

The next section will detail some physical distributions 
for the two models discussed previously. 



(a) 




-1.0 -0.5 0.0 0.5 1.0 

COS0 




COS0 



Figure 8. Angular distributions for fermion and boson pro- 
duction through a resonant graviton. The graviton has a mass 
of 1 TeV. The black line denotes the analytical result and the 
crosses show the simulation data for (a) gg — > G — > e + e~ (b) 
uu — > G — > e + e~ and (c) uu — » G — > 77. 



4 Results 

The following sections will show some distributions pro- 
duced using the new BSM code in Herwig++. As there are 
currently only two models implemented these will form the 
basis of the distributions considered. 



4.1 Graviton Resonances 

The LHC may give us the possibility of detecting nar- 
row graviton resonances at the TeV scale through vari- 
ous hard subprocesses. To test our implementation of the 
RS Model we have picked three processes involving gravi- 
ton exchange, gg — > G — > e + e _ , uu — > G — > e + e~ and 
uu — > G — > 77. The plots of the angular distribution of 

6 This includes all of the spin combinations needed in the 
MSSM and RS model. Additional cases will be implemented 
when required for a new type of model. 



the outgoing fcrmion/boson with respect to the beam axis 
in the centre-of-mass frame are shown in figure [H There is 
good agreement here with the analytical result from the 
matrix element and the numerical simulation indicating 
the correct implementation of the graviton Feynman rules 
and new matrix elements. 

These distributions show the characteristic behaviour 
of an exchanged spin-2 particle. The angular dependence 
of an exchanged spin-1 boson on the other hand is no- 
tably different and therefore this kind of distribution is 
extremely useful in identifying the two cases and eliminat- 
ing possible background spin-1 exchange when searching 
for this new mode [17] in future experiments. This kind of 
behaviour may be important for the LHC since Randall- 
Sundrum type models predict the possibility of narrow 
graviton resonances at the TeV-scale [18]. Discovery of 
the kind of behaviour shown in figure [8] would certainly 
be a very strong indication of the existence of some type 
of extra dimensions model. 
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4.2 Squark Decay 

The spin correlation algorithm discussed in section [2] was 
shown to work in the case of ti production and decay. 
One of the simplest cases to consider for a SUSY model 
is the different decay modes of a left-handed squark. Con- 
sidering the decay of the squark via the two modes (a) 

Ql — > X2I ~ * IrI +u and (b) q_L — ► X%1 — * IrI~ u and 
plotting the mass distribution of the produced quark and 
(anti-)lepton allows the effect of spin correlations to be 
shown. 

The plots in figure[9]were produced at Snowmass point 
5 [19] where tan/3 = 5, sign(^i) = +, M = 150 GeV, 
M 2 = 300 GeV and A = -1000 GeV. This parameter 
set gives the following particle spectrum using SOFTSUSY 
2.0.8 [20] Mu L = 672.82 GeV, = 231.29 GeV and 

Mr = 192.87 GeV. 

There is a stark difference in the quark-lepton mass 
distribution for the two decay modes considered above. 
The difference is due to the helicities of the external par- 
ticles. At the mass scale of the squark the quark can be 
considered massless and left-handed while the produced 
lepton and anti-lepton will be right-handed. When back- 
to-back the lepton-quark system will have net spin-1 and 
as such can not be produced in a scalar decay while the 
anti-lepton-quark system will have spin-0 and is able to 
proceed. 

The end-point in both of the distributions is due to a 
kinematical cut-off where the invariant mass of the quark- 
lepton pair is at a maximum. The value of this end-point 
can be calculated by considering the mass-squared when 
the pair is back-to-back. The value is given in [21] as 

(m? - m? ) (m 2 - mf ) 

(™)L = ^ — ^L- (12) 

m n 

Using the values for the sparticles above one finds a value 
for the cut-off of 348.72 GeV which is consistent with the 
plots in figure El 

4.3 Gaugino Production 

Supersymmetry predicts the existence of Majorana fcrmions 
and it is necessary to ensure that their spin correlations 
are implemented correctly in our new framework. We con- 
sider three production processes and the angular distribu- 
tions of the leptons produced in the subsequent decays. 
The SUSY spectrum for each was again generated using 
SOFTSUSY and the masses for the points used in each 
process are given in the relevant section. 

4.3.1 e+e- -> X2X1 

Here we consider the production of the lightest and next- 
to-lightest neutralinos with the x\ decaying via the two 

modes (i) X °2 - - l+l-xl, (*) xl - - 

I l + Xi at SPS point lb. The relevant sparticle masses are 



M Y o = 306.55 GeV, M v o = 161.78 GeV, Mr = 253.82 GeV. 

A2 1 Xl 1 hi 

Figure [10] shows how the polarisation of the beam affects 
the angular distribution of the lepton produced from the 
X°> decay. 

The lepton shows a correlation with the beam polarisa- 
tion due to the neutralino being a fermion and preserving 
spin information when decaying. Figures [TTJ and [12] show 
the angular dependence of the final-state lepton for the 
case of an intermediate Z° boson and Ir respectively. As 
is to be expected for an intermediate slepton the incoming 
beam polarisation has little effect on the angular distribu- 
tion of the final-state lepton due to its scalar natur^H The 
plots are in good agreement with the HEKWIG results. 

4.3.2 e+e- -> xtxT 

We now consider the production of chargino pairs and 
their associated decays. Two possible decay modes of the 
xt are (a) xt -> W+Xi and (b) xt -> vj^. Here we 
use the W decay mode for the lightest chargino and the 
sneutrino decay mode for the heaviest chargino in or- 
der to consider final-states with differing spins. The mass 
spectrum was generated for SPS point la where = 
377.39 GeV, M % + = 181.53 GeV, M 0l = 185.42 GeV and 
M y o = 97.00 GeV. 

Figure [T3] shows the angle of the produced electron 
for the production of the lightest chargino. As is to be 
expected the beam polarisation affects the lepton distri- 
bution because of the intermediate W boson carrying the 
spin correlations through to the final-state. The effects are 
similar in the case of the sneutrino decay of the heaviest 
chargino shown in figure fMl The lepton accompanied with 
the vl still shows correlations with the beam on account 
of the chargino being a fermion. If the lepton had come 
from the decay of a scalar then there would have been no 
such correlation. 



4.4 Tau Decays 

4.4.1 One Prong Decays 

The tau has a number of leptonic and hadronic decay 
modes. A more detailed analysis of these decays shows in- 
teresting features in the distribution of energy to the decay 
products. A typical tau decay involving several mesons has 
the form t+ — > (nm+)(qm°)v T where nra il denotes n > 1 
charged mesons, i.e. the number of prongs, and qui de- 
notes q > neutral mesons. Here we will consider the 
one prong decay t+ —> p^v T —> n w°v T where the r is 
produced from the decay of a f% . Figure [15] shows our re- 
sults for the fraction of visible energy carried away by the 
charged meson in the two cases where the fx is (a) 100% 
left-handed and (b) 100% right-handed. 

There is a stark difference in the energy distribution 
for two possible mixings of the f 1 in figure [T5] due to the 

7 There is some residual effect due to the correlation of the 
Ir direction with the beam in the X2 decay. 
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(a) (b) (c) 




Figure 10. The angle between the lepton produced in e + e~ — > X2X1 ~~ y Ir}~ an d the beam in the lab frame for a centre-of-mass 
energy of 500 GeV and (a) unpolarised incoming beams, (b) negatively polarised electrons and positively polarised positrons 
and (c) positively polarised electrons and negatively polarised positrons. The black histogram is from HERWIG and the crosses 
from Herwig++. 



resulting helicity of the decaying p. For the case where 
the fi is entirely tl the p has a higher probability of be- 
ing transversely polarised, from the results of [25], which 
favours the equal splitting of energy between the two pi- 
ons as confirmed by the first plot. A t\ that is entirely tr 
however, will give rise to mostly longitudinally polarised 
p mesons that prefer to distribute their energy unequally 
and favour a distribution where one meson receives most 
of the visible energy from the r decay. This is again con- 
firmed in our second plot. The Herwig++ results are plot- 
ted together with those from HERWIG with the TAUOLA 
decay package [24] that is designed specifically for the de- 
cay of polarised r leptons. 

4.4.2 Squark Decay 

The use of the effects described above in the study of 
SUSY models has long been recognised. In [23] a mech- 



anism for the determining the spin properties of parti- 
cles involved in SUSY cascade decays using r polarisa- 
tion was suggested. The method involves analysing invari- 
ant mass distributions of different particle pairs along the 
decay chain q a -> q\l -> t±t^ -> r f T r Xi wit h the r 
decay restricted to — > n v r . The various normalised 
invariant mass distributions are shown in figure [16] for 
q Q = q L at SPS point la where M QL = 558.4 GeV, M x o = 
180.96 GeV, M fl = 134.56 GeV and M x o = 97.00 GeV. 
Since an experiment would be unable to distinguish a near 
or far t/tt their distributions are combined. The normali- 
sation is to the maximum of the invariant mass. For m TT 




which is equal to (m^ 7r ) ma x- In the case of the qr plots the 
maximum of 
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(a) (b) (c) 




0123 0123 0123 

fl e -r fl e -r e e-r 

Figure 11. The angle between the lepton produced in e + e~ — > X2X1 ~~ > Z°Xi ~* l~l + Xi an d the beam in the lab frame for a 
centre-of-mass energy of 500 GeV and (a) unpolarised incoming beams, (b) negatively polarised electrons and positively polarised 
positrons and (c) positively polarised electrons and negatively polarised positrons. The black histogram is from HERWIG and 
the crosses from HerwigH — K 



(a) (b) (c) 




Figure 12. The angle between the lepton produced in e + e~ — > X2X1 ~~ * — * ^~' + Xi an d the beam in the lab frame for a 
centre-of-mass energy of 500 GeV and (a) unpolarised incoming beams, (b) negatively polarised electrons and positively polarised 
positrons and (c) positively polarised electrons and negatively polarised positrons. The black histogram is from HERWIG and 
the crosses from Herwig++. 



(^l-^ieX 1 -^/^)' 
( m L ~ m h) ( 1-m !^ TO ^) ' 

is taken and again this equals (m^) max . 

The differences in shape of the charge conjugate plots 
in figure [H] for the r and ir are due to the different he- 
licities of the r and tt as explained in section 14.21 The 
kinks in these distributions show the change from near 
to far leptons or pions making up the main components 
of the event. These kind of invariant mass distributions 
serve as a good indication of the spin properties of the 
particles involved in cascade decays. This information is 
important when trying to confirm an exact model of new 
physics since it is possible for two different BSM models 
to imitate each other in certain decays even though the 
new particles introduced into each model have different 
spin assignments [26]. 



Again our results are plotted along with those from 
HERWIG using the TAUOLA package. There is excellent 
agreement between the two sets of results and the distri- 
butions follow those of figure 3 in [23] . 

5 Conclusions 

We have described a new method for including new physics 
models in the Herwig++ event generator which is more 
general than the previous approach. It allows new models 
to be implemented with a minimal amount of work. For 
the models implemented the results are in good agree- 
ment with either analytical answers or those from the 
HERWIG event generator. Any new model will automat- 
ically have spin correlations included since the algorithm 
demonstrated in section [2] is formulated independently of 
any specific model and has been shown to agree with ex- 
pected results. 
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(a) (b) (c) 




Figure 13. The angle between the lepton produced in e + e~ — > X1X1 ~ * W + W~XiXi ~ * l~l vi&i and the beam in the lab 
frame for a centre-of-mass energy of 500 GeV and (a) unpolarised incoming beams, (b) negatively polarised electrons and 
positively polarised positrons and (c) positively polarised electrons and negatively polarised positrons. The black histogram is 
from HERWIG and the crosses from Herwig++. 

(a) (b) (c) 



0.6 I — 1 1 1 1 | 1 1 1 1 | 1 1 1 1 _| — 




Figure 14. The angle between the lepton produced in e + e~ — ► X2X2 ~ * and the beam in the lab frame for a 

centre-of-mass energy of 1 TeV and (a) unpolarised incoming beams, (b) negatively polarised electrons and positively polarised 
positrons and (c) positively polarised electrons and negatively polarised positrons. The black histogram is from HERWIG and 
the crosses from Herwig++. 



(a) (b) 




Figure 15. Energy fraction, z, carried away by the charged meson in the one prong r decay p — > ir n° for (a) left-handed fi 
and (b) right-handed fi . The black histogram shows the results from HERWIG with the TAUOLA [24] decay package and the 
crosses indicate the Herwig++ results. 
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In the future we plan to implement other BSM physics 
models in Herwig++ which will enable comparative stud- 
ies of the phenomenological consequences of these models 
to be carried out under the framework of the same gener- 
ator. The current release, 2.0, of HerwigH — h [2] does not 
include any of the features discussed in this paper, it will 
however appear in version 2.1 of the event generator. 

We would like to thank our callaborators on the Herwig+H- 
project for many useful comments. This work was supported 
in part by the Science and Technology Facilities Council and 
the European Union Marie Curie Research Training Network 
MCnet under contract MRTN-CT-2006-035606. 



A Helicity Code 

In the FORTRAN HERWIG program the helicity ampli- 
tudes were calculated using the formalism of Kleiss and 
Stirling [27]. However this meant that: 

— all the matrix elements had to be calculated in the 
same frame; 

— it was impossible to interface the correlations in the 
parton shower and those in the decay of unstable fun- 
damental particles; 

— each matrix element had to be separately calculated. 



In HerwigH — h we choose to use to use an approach based 
on the HELAS formalism for the calculation of matrix el- 
ements. This approach has a number of advantages: 

— we can use the spinors and polarisation vectors calcu- 
lated when the particle is produced to calculate their 
decays in a different frame, after an appropriate Lorcntz 
transformation, so that each step of the calculation can 
be done in the most relevant frame; 

— more complicated matrix elements can be calculated 
from the basic building blocks rather than coded from 
scratch; 

— the inclusion of particles other than scalars, spin-i 
fcrmions and massless spin-1 bosons, which is compli- 
cated in the Kleiss and Stirling formalism, is relatively 
simple. 

The implementation of the HELAS formalism in Herwig+H- 
is based on two fundamental types of objects, WaveFunc- 
tionBase and VertexBase. 

The WaveFu notions contain the momentum and a pointer 
to the properties of the particles together with the basis 
state for a given particle spin.WaveFu notion Base stores the 
momentum of the particle and a pointer to the Particle- 
Data object for the particle. The inheriting ScalarWave- 
Function, SpinorWaveFunction, SpinorBarWaveFunction, Vec- 
torWaveFu notion, RSSpinorWaveFunction, RSSpinorBarWave- 
Function and TensorWaveFunotion classes then contain stor- 
age of the wavefunctions for spin-0, -i, -1, -| and -2 parti- 
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cles together with methods to calculate the wavefunctions 
for a given helicity state. 

The Vertices contain methods to combine the wave- 
functions for a specific vertex to give either off-shell wave- 
functions, which can be used in further calculations, or the 
matrix element. The VertexBase class contains storage of 
the particles which interact at a given vertex. A number of 
classes then inherit from this class and implement the cal- 
culation of the matrix element and off-shell wavefunctions 
for a given Lorentz structure of the vertex in terms of ar- 
bitrary couplings which are calculated by virtual member 
functions. 

This strategy is essentially identical to that adopted in 
the original HELAS approach. However, the new structure 
has the benefit that the user need only decide which ver- 
tex to use with the structure then supplying the relevant 
couplings, whereas these had to be specified by the hand 
when using the HELAS library. 



A.l Conventions 

To numerically evaluate the matrix elements using the 
HELAS formalism we need a specific choice of the Dirac 
matrices, we currently support two options. The conven- 
tional low-energy choice, used in for example [28] 



HABER _ / °i 

1 -<Ji o 



HABER 
7o 



1 

-1 



HABER _ / 1 
'5 — 1 1 



(13) 



and the original choice of HELAS which is more appropri- 
ate at high energies 



7< 



HELAS _ 



a, 

-<Tt 



.7o' 



HELAS 



1 

1 



75 



HELAS _ 



-1 
1 



(14) 



These two representations are related by the transforma- 
tion 

^helas = S^Habcr where 5=-^=^j j 1 ^ . (15) 

A number of container classes are implemented in the 
ThePEG framework [29], on which Herwig++ is built, to 
store the basis states for the different spins: 

— LorentzSpinor: storage of the spinor, u or v, for a spin-i 
fcrmion; 

— LorentzSpinorBar: storage of the barred spinor, u or v, 
for a spin-i fermion; 

— LorentzPolarizationVector: storage of the polarization 
vector, e^ 1 , for a spin-1 boson; 

— LorentzRSSpinor: storage of the spinor, u M or v^, for a 
spin- 1 fermion; 

— LorentzRSSpinorBar: storage of the barred spinor, u M 
or w M , for a spin-| fermion; 



— LorentzTensor: storage of the polarization tensor, e^ u , 
for a spin-2 boson. 

In addition to providing storage of the basis state these 
classes implement the Lorentz transformations, both boosts 
and rotations, for the objects. In the case of fcrmions the 
Dirac basis used to calculate the spinor, together with 
whether the spinor is u or v type is also stored. Methods 
to convert between the two supported Dirac matrix defi- 
nitions are implemented together with the transformation 
between u and v spinors for Majorana particles. 



A. 2 Lorentz Transformations 

In addition to the storage of the basis states we need to be 
able to transform them between different Lorentz frames. 
The Lorentz transformation for a spinor is given by 



^\ x ') = tf{ax) = S(a)ip(x), 



(16) 



where 



dx'" 



For a Lorentz boost along the direction 



specified by the unit vector n with a magnitude [5 the 
transformation is given by 

S boost = cosh (I) + sinh (|) n l7 V , (17) 

where tanhx = P- For a rotation by an angle <j> about the 
unit vector n the Lorentz transformation is given by 



^rotation — COS 



e^n fc7 y. (18) 



2 ) \2 

The Lorentz transformations for a four vector is given by 

e*V) = L{ay v e"{x). (19) 

If we wish to boost by a factor (3 along a unit vector n the 
transformation is 



Lt = 



7 -7/3ni 
— 7/3ni 1 — fiihiuj 
-7/3n2 — nih\ij) 1 
-7/3n 3 - h^fiiuj 



-Jf3fl2 - 

- h\Jl2LO — 

- h^hiu — fi2h 3 w I 

- n 3 n 2 uJ 1 - n 3 n 3 uj ) 



j(3n 3 \ 
™ W I , (20) 



where ui = 1 — 7 and 7 



■■. The Lorentz transfor- 



mation for a rotation by an angle S about a unit vector n 
is 

Lt = 

/l 

\h\fii + c§ \fi1fi2 — s$n 3 \fi1h3 + sgfi2 



An 2 ni + sgn 3 An 2 n 2 + c$ \h2fi3 — sghi 
\0 Xh 3 ni — ssh 2 Xn 3 h 2 + s$hi Xh 3 h 3 + c$ 



(21) 



where cs = cos 5, s$ = sin<5 and A = 1 — cos 8. The trans- 
formations for the higher-spin particles can then be con- 
structed as products of the spin-i and spin-1 transforma- 
tions, i.e. 



r(x') = L(a)ZS(aW(x), 



e^(x') = L{aY a L{aT p e^. 



(22a) 
(22b) 
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A. 3 WaveFunctions 

A. 3.1 ScalarWaveFunction 

The ScalarWaveFunction class inherits from the WaveFunc- 
tionBase class and implements the storage of the wavefunc- 
tion of a scalar particle as a complex number. For exter- 
nal particles this is just 1, however it can assume different 
values when the WaveFunction is the result is an off-shell 
internal line from a Vertex class. 



A. 3. 2 SpinorWaveFunction and SpinorBarWaveFunction 



where m is the mass and p T = ^jp\ + p y . We include two 
choices of the polarization vectors 



e"(p,A = ±1) 



1 



V2 

e"(p,A = 0) = e£(p) 



(TeSM-ieSfo)), 



(27a) 
(27b) 



which is the choice used in HELAS. However, while this 
option is available in the Herwig++ by default we include 
the additional phase factor exp(iA0) as in [28] in order 
to make the inclusion of spin correlations in the parton 
shower easier. 



As with the ScalarWaveFunction the SpinorWaveFunction 
and SpinorBarWaveFunction classes inherit from the Wave- 

FunctionBase class. The spinor is stored as either a LorentzSpirA>.r3.4 RSSpinorWaveFunction and 
or LorentzSpinorBar. In addition the calculation of the spinors RSSpinorBarWaveFunction 

for external particles is implemented. The spinors are cal- 
culated in terms of two-component spinors, as in [16] 



x+{p) 



x-(p) 



1 



v^IpKIpI 
i 



+ Pz) 



V 2 \p\ (IpI +p* 



Ip 

Px - 

-Pa 
|pj 



+ Pz 
+ iPy 

+ iPy 
+ Pz 



(23a) 



(23b) 



where p x , y ,z are the x, y and z components of the momen- 
tum respectively, E is the energy of the particle and |p| 
is the magnitude of the three momentum. For the HELAS 
choice of the Dirac representation the spinors are given by 



u(p) 



v(p) = 



u-\{p)x\{p) 

u\(p)x\(p) 

^\{p)x-\{p) 
\(p)X-\{p) 



where u>±(p) — \f~E ± |p| and the helicity A 
larly for the low energy definition 



= ±1. 



u{p) = 



v(p) 



^+(p)x\(p) 
^-(p)x\(p) 

u-(p)x-\(p) 
-\uj+(p)x-\(p) 



Although there are currently no fundamental spin-| par- 
ticles included in Herwig++ the Rarita-Schwinger spinors 
for spin- 1 particles are included both to allow the sim- 
ulation of spin- 1 hadronic resonances and for the possi- 
ble future inclusion of the gravitino. The RSSpinorWave- 
Function and RSSpinorBarWaveFunction inherit from the 
WaveFunction Base class and implement the storage of the 
Rarita-Schwinger spinors using the LorentzRSSpinor and 
LorentzRSSpinorBar classes respectively. 

The spinors are calculated using the Clebsch-Gordon 
decomposition: 



(24a) 




-2) = 


e"(p,-lty(p,-l); 


(28a) 


(24b) 


V{pA = 


-1) = 


yle"(p,-lMP,l) 


(28b) 


Simi- 






+ y|^(p,0)^(p,-l); 


(28c) 


(25a) 




= 1) = 


y|^(p,l)V(p,-l) 


(28d) 


(25b) 


r(p,x 


= 2) = 


+^(P,0); 
e"(p,l)^(p,l). 


(28e) 
(28f) 



A. 3. 3 VectorWaveFunction 

The VectorWaveFunction class inherits from the WaveFunc- 
tion Base class and implements the storage of the polariza- 
tion vector using the LorentzPolarizationVector class. The 
polarization vectors of a spin-1 particle can be calculated 
using, as in [16] 



For massive particles the spinors are calculated in the rest 
frame of the the particle and then boosted to the required 
frame in order that the Clebsch-Gordon decomposition 
can be easily applied. For massless spin-| particles, which 
only have the ±2 helicity states, the spinors are calculated 
in the same frame as the momentum. 



£i(p) = T~\ — ( Q >PxPz,PyPz, -Pt) 
\P\Pt 



e 2(p) = {0,-Py,Px,0), 

PT 

e 3(W = [ —iPx^PyiP 



m|p 



(26a) 
(26b) 
(26c) 



A. 3. 5 TensorWaveFunction 



The TensorWaveFunction class inherits from the Wave- 
Function Base class and implements the storage of the po- 
larization tensor for spin-2 particles using the LorentzTen- 
sor class. 
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The wavefunction is calculated using the Clebsch- 
decomposition: 



A = 


-2) = 


e"(p,-l)e"(p,-l); 


(29a) 


e^iv A = 


-1) = 


jl \ e n( p -l) e v {v 0) 


(29b) 






+ e"(p,0)e"(p,-l)]; 


(29c) 




= 0) = 


y|[ e ^,l)e>,-l) + e <>, 


-l)c"(p,l) 






+2e"(p,0)e v (p,0)]; 


(29d) 


e^(p,A 


= 1) = 


y|[^(p,l)e"(p,0) 


(29e) 






+^( P ,0)e"( P ,l)]; 


(29f) 




= 2) = 


£*>,!)£>, 1). 


(29g) 



Here this is applied in the frame in which the momentum 
is specified. 

A. 4 Vertices 

The Vertices all inherit from the VertexBase class. In gen- 
eral for all the vertices all the particles and momenta are 
defined to be incoming. 

A. 4.1 Scalar Vertices 

There are a number of vertices involving scalar bosons. 

FFSVertex The vertex for the coupling of a fermion and 
antifermion to a scalar boson is defined to have the 
perturbative form 

*c/ 2 a A P A /i03, (30) 

where c is the overall normalisation, a A are the left /right 
couplings, P\ are the helicity projection operators, fx 
is the wavefunction for the fermion, f 2 is the wavefunc- 
tion for the antifermion and tp 3 is the wavefunction for 
the scalar boson. 
GeneralSVVVertex In addition to the perturbative form 
for the vertex coupling a scalar and two vector bosons 
described below we include a general form for this in- 
teraction so that effective vertices, for example h° — > 
gg, can be implemented. The form of the vertex is 

ic [aoog^ + a 22 p^ + d2iP 2 P 3 + a 3 2P 3 P 2 + a 33 P 3 P 3 
+a e e^p 2a p 3/3 ] e 2 „e 3v <t> 1 , (31) 

where p2,3 are the momenta of the vector bosons, 62,3 
are the wavefunctions of the vector bosons, <f>\ is the 
wavefunction of the scalar boson, c is the overall cou- 
pling and dij are the couplings of the different terms. 
SSSVcrtex The vertex for the coupling of three scalar 
bosons is defined to have the perturbative form 

icfafcfa, (32) 

where ^1,2,3 are the wavefunctions for the scalar bosons 
and c is the coupling. 



SSSSVertex The vertex for the coupling of four scalar 
bosons is defined to have the perturbative form 

««/>i</>20304, (33) 

where (^1.2,3.4 are the wavefunctions for the scalar bosons 
and c is the coupling. 
VSSVertcx The vertex for the coupling of a vector boson 
and two scalar bosons is defined to have the perturba- 
tive form 

- ic(p 2 -p 3 ) ■ ei02^3, (34) 

where t\ is the wavefunction of the vector boson, </> 2 ,3 
are the wavefunctions for the scalar bosons and j>2,3 are 
the momenta of the scalar bosons and c is the coupling. 
VVSSVertex The vertex for the interaction of two vector 
and two scalar bosons is defined to have the perturba- 
tive form 

icg^ei^e^fofa, (35) 

where ei )2 are the wavefunctions of the vector bosons 
and 3i 4 are the wavefunctions for the scalar bosons 
and c is the coupling. 
VVSVertex The vertex for the interaction of two vector 
and a scalar boson is defined to have the perturbative 
form 

icg^e^e^fo, (36) 

where ei ;2 are the wavefunctions of the vector bosons 
and </>3 is the wavefunction for the scalar boson and c 
is the coupling. 

A. 4. 2 Vector Vertices 

There are a number of vertices involving vector bosons. 

FFVVertex The interaction of a fermion, antifermion and 
a vector boson is taken to have the perturbative form 

icf 2 ^a x P x he 3lM (37) 

where c is the overall normalisation, a A are the left /right 
couplings, fi is the wavefunction for the fermion, f 2 
is the wavefunction for the antifermion and £3 is the 
wavefunction for the vector boson. 
VVVVertex The interaction of three vector bosons is taken 
to have the perturbative form 

ig[(pi -P2Vg a0 + ip2-p 3 ) a g^ 

+ (P3-Pi)V 7 ]eiae2/3e 37 , (38) 

where £1,2,3 are the wavefunctions of the vector bosons 
and pi,2,3 are the momenta of the vector bosons. 
VVVVVertex The interaction of four vector bosons is taken 
to have the form 

ic 2 [2ei • £ 2 £ 3 • £4 - £1 • £3£2 • £4 - £1 • £4£2 ■ £3] , (39) 

where £1,2,3,4 are the wavefunctions of the vector bosons 
and Pi, 2,3,4 are the momenta of the vector bosons. For 
the quartic gluon vertex this is the contribution of one 
colour structure. The others can be obtained by an 
appropriate reordering of the input wavefunctions. 
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A. 4. 3 Tensor Vertices 

There are a number of vertices involving spin-2 particles. 
The form of the Feynman rules follows that of [22]. 

FFTVertex The interaction of a pair of fermions with a 
tensor is taken to have the perturbative form 



- f h [7m (Pi - Vi)u + IviPl - P2)p 
-2flW(yi ^ #2) + Ag^rrif] fie% u 



(40) 



where re is the defined as 2/ A cut - f /, pi t 2 are the mo- 
menta of the fermions, f\ is the fermion wavefunction, 
fi is the antifermion wavefunction and eg" is the po- 
larisation tensor for the spin-2 particle. 
VVVTVertex The interaction of three vector bosons with 
a tensor is taken to have the perturbative form 

fff [Cpu,pa{Pl ~P2)\ + C^.pxipz -p\) a 

+ C^ a x {P2-P3)p + F MV ,pax] 4^44" (41) 

where re is 2/ A cut _ ff, are the momenta of the 

vector bosons, e± 2 3 are the polarisation vectors and 
e^ v is the polarisation tensor. The C and F symbols 



are defined as 



a 



p,v,pa 



9p,p9va T QpaQvp QpvQpm 



Fp U , pa \ = g w g a \{p2 ~ p-i)v + gpag P \(p3 ~ Pi) 
+gp\g P a(pi - P2)u + 0« <->• v). 



(42) 



(43) 



VVTVertex The interaction of two vector bosons with a 
tensor is taken to have the perturbative form 



- [{ml + p x ■ p 2 )C fJtv , ptx + D„] cfef ef , (44) 



where re is the defined as 2/ A cu t- ff, mv is the mass 
of the gauge boson pi.2 are the momenta of the vector 



bosons, e{ 



2..S 



are the polarisation vectors and e 3 



' 1 



is the 



polarisation tensor. The C symbol is defined as above 
and D is defined as 

Dp^pa = gpvkiak 2 p - [gpaki v k 2p + gppki a k 2u 

-g pa k\pk 2v + (n <-> v)} . (45) 

FFVTVertex The interaction of a pair of fermions with a 
vector boson and a tensor is taken to have the pertur- 
bative form 



ign 



fiiC^pa - gp V g P a)l a h4 e T 



(46) 



where re is the defined as 2/ A cut - f f, £3 is the polar- 
isation vector for the boson, j\ is the fermion wave- 
function, f 2 is the antifermion wavefunction and e^" 
the polarisation tensor for the spin-2 particle. The C 
symbol is defined above. 
SSTVertex The interaction of a pair of scalars with a ten- 
sor is taken to have the perturbative form 



IK 



[m 2 s g pv - p\pPi v ~ P\ v P2p + g^vPi ■ P2] ^ v (f>i4> 2 

(47) 



Class Name 


Process Type 


MEff2ff 


T t' t" t'" 


MEfT2ss 


<P <P — » <f> <f> 


MEff2vv 


& <P -> V V 


MEfv2fs 


$ V -> W </> 


MEvv2ss 


v v -> <j> </> 


MEvv2ff 


V V' w 


MEfv2vf 




MEvv2vv 


1 n in 

VV V V 



Table 1. The MatrixElement classes. 



Class Name 


Decay Type 


FFSDecayer 


& — ► w cj> 


FFVDecayer 


V^rW' V 


SFFDecayer 


<f> — > >i> i> 


SVVDecayer 


<i>-> V V 


SSSDecayer 


1 11 

</>—><!> <f> 


SSVDecayer 


(j>-xj>V 


TFFDecayer 


T -> & W 


TSSDecayer 




TVVDecayer 


T^VV 


VFFDecayer 


V -» $ 9 


VSSDecayer 




VVVDecayer 


1 II 

V -> V V 



Table 2. The two-body Decayer classes. 



where re is the defined as 2/ A cut - f f, A cut - ff is the 
ultraviolet cutoff scale, ms is the mass of the scalar, 
Pi ; 2 are the momenta of the scalars, the polari- 
sation tensor for the spin-2 particle and (j>i t 2 are the 
scalar wavefunctions. 



B MatrixElement Classes 

Table [1] gives a list of the implemented 2^2 Matrix- 
Element classes with a description of the types of hard 
subprocess they can handle. \P stands for a fermion, V a 
vector boson and cj> a scalar. Charge conjugate modes are 
handled by the class responsible for standard mode. 



C Decayer Classes 

Table [2] is a list of the implemented two-body Decayer 
classes with a description of the types of decay that they 
are designed to handle. W stands for a fermion, V a vector 
boson, <j) a scalar and T a tensor. Charge conjugate modes 
are handled by the class responsible for standard decay 
mode. 
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icr" ( Pl +p 2 T 



Figure 17. Feynman rule for the interaction of a gauge boson 
with a pair of sfermions. The definition of P for the various 
types of gauge boson and sfermion is given in table [3] The 
momenta are to be taken as in the direction of the arrows. 



D Implemented Vertices 

Below is a list of vertices and associated Feynman rules, 
for BSM physics, as they are currently implemented. The 
rules involving colour are written with the colour depen- 
dence explicitly pulled to the front since this is not in- 
cluded in the code for a specific vertex as explained in 
section [231 If a structure constant f abc is involved an ad- 
ditional factor of i is pulled out along with it due to the 
commutator relation [i a ,i 6 ] = if a 



pabc-^c 



D.l RSModel 

The Feynman rules as implemented in HerwigH — h are given 
in section [A. 4. 31 



D.2 MSSM 

We give the Feynman rules for the MSSM as implemented 
in HerwigH — \-. The sfermion mixing matrices are denoted 
by Q£ and L k aQ for the squarks and leptons respectively 



J a{3 

where k is the generation number, a the left /right eigen- 
state and (3 the mass eigenstate. 2Vy, Uij and Vij are the 
neutralino and chargino mixing matrices respectively. The 
primed matrices in the neutralino rules are related to the 
unprimed ones via 



Nn cos Qw + N i2 sin 9 W , (48a) 

N i2 cos W - Nu sin 9 W , (48b) 

N l3 , (48c) 

N i4 . (48d) 



z° 



z° 



z° 



z° 



\y- 



ff 



qq 
If 



qq 



U a U0 



d a d 



ad/3 



la I 



q a qp 



vlr- 



J. ( I 



cos 6\y 



-g 

_J2_ 



-B 



pff 

1 a/3 

7 — 

sqq 

°a(3 
T^T — 

a' ' 

°af3 



s qq 

°af3 



\{-Qf a Q% + 2e u sm 2 e W 5 a g) 



I {Ll^ 1 ^ 1 -2sm 2 6 w S a0 



■hSn 



Q2i /-\'2i—\ 



L 



1/3 



Table 3. Couplings for the gauge bosons and sfermions. 

•/i(pi) 

ic7 M [a L P L + a R P R ] 

fj(P2) 

Figure 18. Feynman rule for the interaction of a gauge boson 
and a pair of gauginos. The couplings are defined in table [4] 



V 




X?, 9 




ic [o, l Pl + a R P R ] 



Figure 19. Feynman rule of the coupling of the neutralinos 
and gluino with a Standard Model fermion and an sfermion. 
Table [5] gives the definitions of ol and a R . 



fi fia 


c 


a L 


a R 


Ui did 


-g 


m «i v /2 r>2i-i 


rr r^2i-\ ' n d i °'j2 ^,2i-\ 
^Jl^lc v'lMw cos V 2c 




-g 




T T r<2i rn Ui V j2 ^2i 


^ilia 


-g 





rr r2i-l m ( i Lr j2 r 2i-l 


f>i Via. 


-g 


m h u h 

V2M lv cos 





Table 6. Chargino couplings for the i th fermion generation. 
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VP 
V 


J J 






ft 


7 


Xi Xj 


— e 


1 


1 


W+ 


m 


9 






z° 


Xi Xj 


9 

cos 9 W 


— \Na N* 3 + \ Nu Nj 4 


-al 


z° 


Xi Xj 


9 

COS By/ 


-VilV*! - \Vi2V* 2 

+Sij sin 2 6w 


—Uullji — ^U* 2 Uj2 
+Sij sin 2 9w 


9a 


9b 9c 


{ifabc} g s 


1 


1 



Table 4. Feynman rules for the coupling of a gauge boson to a pair of electroweak gauginos. All momenta are to be taken as 
outgoing. 



fi }" 



d% di, 



li li. 



V2 



-V2 



-V2 



-V2 



{t a bc } (-g.y/2) 



a L : 



g m "i Af j4 n 2i 
2M W sin /3 Via 



cos 0\y 



2M W sin / 



Qf a + via ee 



M l . 9l|-e» sin 2 8 w )N' j2 \ 
" JV jl "r cos% I 



9 m dj N j3 rfii-l _ n 2 ^ 1 ( ee AT'* — 9 " d S '"^ flw,JV j2 
ai -- 2M W cos/3 Via V2a I eedlVji cos W 



Or 



. 9m di N j3 n 2i-l , n 2i-l ( N l _ 9(^+e d sin 2 9 w )N' j2 \ 
■ 2M w cos/3V2a ' Via I C W jl cos 6 W J 



. g" 1 !,^ r2»-l , r2»-l / v /* ffsi 
ttt ' 2MhTcos^ la + ^2a [ e ^jl ~ ~ 



2M„,cos/3 Jj 2a 



- ^1 



a L : 



9N' j2 
an ■ 2 cos e n 



1 (eNj 



Sin gjy jV" 7 * 

cos 

g (l-sin 2 %) 
Jl "T" cos« w iV j'2 



ax, : 



-(?2a 



a K : Qla 



Table 5. Neutralino and gluino couplings. In the case of the gluino k = 2i for up-type quarks and k — 2i — 1 for down-type 
quarks. 



fia fi0 


r 


- ~ * 
Uia V>i(3 


Mz c rt +W [Qf a Q%{\ - sin 2 W ) + e u sin 2 6 w Qf a Q%] 

2 

fn u . cos a r 2i 2i 2i 2i 1 
M w sin/3 [QlaQlp + QlaQip] 

-aCrinfl^ cos a + M sin a) [QlaQiJg + 


dia d^p 


M Z sin(a + /3) 

cos 6>v^ 
m 2 , sin a [" 
+ M w cos [Ql 

i md * r \ 

~ 2M W cos /3 y 


Ql^QVi 1 ^ + ed sin 2 ^) - e d sin 2 ^Q^-^^" 1 ] 

a Vl/3 T V2a l <2/3 

i N r/n2i-l/n2i-l . s\2i— l^-fii— ll 

J( sin a + ^ cos a) Q 2 a Q1/3 + Qia Q213 \ 


T. i* 


M z sin(a+/3) f r 2i- 1 r 2i- 1 1 1 „;„2 fl , „;„2 a r 2i-l r 2i-l] 
cos% [ L la ^1/3 U Sln fwj+Sm ^^2a ^2/3 

, ™?iSin" \ T 2i-\ T 2i-\ , r 2i-l r 2i-ll 
"T" M W cos /3 L^la ^1/3 ^2a ^2/3 

+ 2M^(^. sin « + MCosa) [Llj.-^^+L"- 1 ^- 1 ] 


Via Vifi 


M Z sin(a + ,3) 
2 cos 



Table 7. ft couplings to sfermion pairs with c = ig. 
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fitx fi@ 


r 


11. ?/, ■ a 


~ Mz ::: ( el +m [Q 2 ^Q 2 M\ - e « ^ Ow) + e u sin 2 0wQf a Q%] 

M w sinf3 [QlctQip + QlaQw] 

- 2M Z U sinp sin ct cos ct) [Qf a Q% + Qf a Q%} 




Af z cos(a+/3) [^21-1^,21-1/1 , „ „;„2 a \ „ „;„2 a rfli-l 

cos9 w [Wia Qip ( 2 +e d sm 0wJ e d sm 0wy 2a 

m d . cos a 1^21-1^21-1 , s~,2i-l f~,2i-l \ 
M^cos/3 [Via ^1/3 + <^2a <^2/3 I 
i md i 1 ■ A \ \rfli—ls\2i—\ | ^21 — 1^21—1 

+ 2M„,cos/3(A tsma A di cosa) Q 2a Q 1(3 +Q la Q 20 


^2/3 ] 


7 r* 


M z cos( Q+ /3) r L «-l L «-l ( l gin 
COS t?JV loj 1/3 V 2 
2 r 

m i, cosa Tr2i-ir2 l -i , r 2*-l r 

M w cos/3 L^la ^1/3 +-^2a i 
+ 2A/ w cos/3^ Sma ^COSa) 


2 dw) + sin 2 Sivijia J -^2^ 1 

2i-l] 
"2/3 J 

r 2i-l 7-2i-l | r 2i-l t- 2i-l] 
^2a ^1/3 + L lc ^2/3 




Via V*g 


Ai" z cos(a+/3) 
2 cos 0^ 



Table 8. // couplings to sfermion pairs with c = 



■ Ui, Vi 



fii 



ic [cllPl + a R P R ] 



h°,H°,A° 



cF 



v 



i/3 




ic [cilPl + a R P R ] C 



Figure 20. Feynman rule of the coupling of a chargino with 
a Standard Model fermion and an sfermion. C is the charge 
conjugation matrix. Table Ogives the definitions of a_L and a R . 




i{\8ab + d abc t c ) gl 



Figure 21. Feynman rule for the coupling of two gluons to a 
pair of squarks. 



Figure 22. Feynman rule for the coupling of a higgs with 
a pair of sfermions. The definition of P for the three neutral 
Higgs' are given in table [7] HO and [9] 



/iQ fi{3 


r 


~ ~ * 


2M„ ( A Vi COt I 3 + M) <W/3 


dia. dig 


2M ; (M tan/3 + fi)5 a ^i3 


ha lip 


2Mw (^ e - + 



Table 9. A couplings to sfermion pairs with c = g. 
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